Data explorations
First thing we did was read in the main data file and merge it with the day 0 data. Because the day zero are true for all temperatures, we duplicate the data for the 19 and 55 degree conditions. Then, we convert log-scale the the indicator bacteria counts, adding a dummy value of 1 to prevent errors of logging 0.
main_data <- read.csv("./data/clean/trial_data.csv")
day0 <- read.csv("./data/clean/day0.csv", skip=3, skipNul = T, stringsAsFactors = F)
day0 <- day0 %>% select(-fs_conc, -fs_vol, -inoc_conc, -inoc_vol) %>% spread(key = measurement, value = Value)
day0<- rbind(
day0,
day0 %>% transform(Temp=replace(Temp,Temp==37, 19)),
day0 %>% transform(Temp=replace(Temp,Temp==37, 55))
)
names(main_data)
## [1] "Temp" "OL" "Time" "Recipe" "Ammonia"
## [6] "pH" "sCOD" "TS" "VS" "Coliforms"
## [11] "Ecoli" "Enterococci" "Methane"
names(day0)
## [1] "Recipe" "Temp" "OL" "Time" "Ammonia"
## [6] "Coliforms" "Ecoli" "Enterococci" "Methane" "pH"
## [11] "sCOD" "TS" "VS"
sn_raw <- rbind(main_data, day0)
sn <- sn_raw %>%
rowwise() %>%
transform(
Recipe=as.factor(Recipe),
Ecoli = log10(Ecoli + 1),
Coliforms = log10(Coliforms + 1),
Enterococci = log10(Enterococci + 1))
summary(sn)
## Temp OL Time Recipe Ammonia
## Min. :19 Min. :0.5 Min. :0.000 FWS12:69 Min. : 865
## 1st Qu.:19 1st Qu.:0.5 1st Qu.:3.000 FWS13:69 1st Qu.:1452
## Median :37 Median :2.0 Median :6.000 FWS31:69 Median :1625
## Mean :37 Mean :2.0 Mean :5.217 Mean :1619
## 3rd Qu.:55 3rd Qu.:3.5 3rd Qu.:6.000 3rd Qu.:1768
## Max. :55 Max. :3.5 Max. :9.000 Max. :2305
## pH sCOD TS VS
## Min. :6.840 Min. : 6522 Min. :4.840 Min. :3.120
## 1st Qu.:7.240 1st Qu.: 8262 1st Qu.:5.570 1st Qu.:3.410
## Median :7.440 Median :10111 Median :5.740 Median :3.560
## Mean :7.436 Mean :10502 Mean :5.874 Mean :3.705
## 3rd Qu.:7.510 3rd Qu.:12261 3rd Qu.:6.035 3rd Qu.:3.860
## Max. :8.060 Max. :22965 Max. :8.095 Max. :5.371
## Coliforms Ecoli Enterococci Methane
## Min. :0.000 Min. :0.000 Min. :0.000 Min. : 0.00
## 1st Qu.:2.249 1st Qu.:2.154 1st Qu.:4.511 1st Qu.: 5.75
## Median :4.772 Median :4.502 Median :5.178 Median : 19.50
## Mean :3.857 Mean :3.675 Mean :4.757 Mean : 40.89
## 3rd Qu.:5.380 3rd Qu.:5.255 3rd Qu.:5.871 3rd Qu.: 75.55
## Max. :8.020 Max. :8.020 Max. :7.322 Max. :188.30
It may prove useful to note the components of the three recipes. The organic loading rates were calcualted based on the volitile solid quantity. the recepies are ratios of 2 of three ingrrerdiates (not including the innoculum, which is added as discussed above): dairy cattle slurry(DCS), Fats Oils and Grease (FOG), and restauraant food waste (RFW). - FWS12: 1:2 FOG:DCS - FWS13: DCS:RFW 1:3 - FWS31: DCS:RFW 3:1
We add those proportion to the data. Then, we save the cleaned, prepated data.
recipes <- data.frame(
Recipe= c("FWS12","FWS13","FWS31"),
RFW= as.factor(round(c(0, 3/4, 1/4), 2)),
DCS=as.factor(round(c(1/3, 1/4, 3/4),2)),
FOG=as.factor(round(c(2/3, 0, 0), 2))
)
sn <- merge(sn, recipes, by="Recipe")
write.csv(sn, "data/clean/cleaned_combined_minitrial_data.csv")
Note that the conditions were set to have 3 different loading rates, which we will treat as categorical because while they are numerical values of .5, 2, or 3.5, these are not observed values. There for we set them to be factors (which allows for categories to be ranked by a inheirent value).
Next, we can do some explortatory analysis:
sn$OL <- as.factor(sn$OL)
ggscatmat(sn, color="Recipe", alpha = .3) + theme(legend.position = "bottom")
## Warning in ggscatmat(sn, color = "Recipe", alpha = 0.3): Factor variables
## are omitted in plot

#ggpairs(sn, aes(color=Recipe, alpha=.5))
Initial observations: - the original Box Behnken design can be seen looking at the top left area, seeing the distribution of values between Time, OL, and Temp. - E. coli and Coliforms are, unsurprisingly, very highly correlated. - TS and VS are highly correlated
Lets replot without some of the highly correlated ones that
ggscatmat(sn %>% select(-TS, -Coliforms), color="Recipe", alpha = .3) + theme(legend.position = "bottom")
## Warning in ggscatmat(sn %>% select(-TS, -Coliforms), color = "Recipe",
## alpha = 0.3): Factor variables are omitted in plot

Lets look at the indocators ### Enterococci
ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) +
geom_point() + facet_wrap(~OL+Recipe) + geom_smooth(method = lm) +
labs(caption=fn("Effect of OL and Recipe"), color="")

It looks like for FWS12 the levels of enterococci stay the same, but in all the FSW13, FWS31, the levels seems to increase or decrease over time. We can try to determine if that is due to
ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) +
geom_point() +
labs(caption=fn("Effect of OL vs DCS"), color="") +
facet_wrap(~OL+DCS) +
geom_smooth(method = lm)

ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) +
labs(caption=fn("Effect of OL vs RFW"), color="") +
geom_point() + facet_wrap(~OL+RFW) +
geom_smooth(method = lm)

ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) +
geom_point() +
labs(caption=fn("Effect of Temp"), color="") +
facet_wrap(~RFW) + geom_smooth(method = lm)

ggplot(sn, aes(x=Temp, y=Enterococci, color=as.factor(RFW))) +
labs(caption=fn("Effect of RFW"), color="") +
geom_point() + geom_smooth(method = lm)
Those rates seem to indicate 2 things:
- that the higher proportion of restaurant food waste, the better enterococci survive
- regardless of the proportion of RFW, increased temperature inhibits Enterococci
E. coli
ggplot(sn, aes(x=Time, y=Ecoli, color=as.factor(Temp))) +
geom_point() +
facet_wrap(~OL+Recipe) +
geom_smooth(method = lm) +
labs(caption=fn("Effect of OL and Recipe"), color="")

E coli counts seem to depend largely on Temp rather than OL or Recipe.
ggplot(sn, aes(x=Time, y=Ecoli, color=as.factor(Temp))) +
geom_point() + geom_smooth(method = lm) +
labs(caption=fn("Effect of Temperature"), color="")
This may indicate that while die-off occurs at cold temperature, higher temperatue increases the die-off rate.
Methane
Methane is the only thing leaving the bioreactors (as biogas), and therefore was measured daily until detructive sampling at 3, 6, or 9 days. We have aggregated the data in the original excel sheet.
ggplot(sn, aes(x=Time, y=Methane, color=as.factor(Recipe))) +
geom_point() + geom_smooth(method = lm) +
facet_wrap(~Temp ) +
labs(caption=fn("Effect of Temperature and Recipe"), color="")

ggplot(sn, aes(x=Time, y=Methane, color=as.factor(Temp))) +
geom_point() + geom_smooth(method = lm) +
facet_wrap(~Recipe ) +
labs(caption=fn("Effect of Temperature and Recipe"), color="")

ggplot(sn, aes(x=Time, y=Methane, color=as.factor(Temp))) +
geom_point() + geom_smooth(method = lm) +
facet_wrap(~Recipe + OL ) +
labs(caption=fn("Effect of Temperature and OL"), color="")
It looks like the rate of poduction is highest at 37, and that at the optimal feed stocks with a high OL yield high prroduction, expecially at 37.
Modelling
stan_data = list(
nbugs=3,
nrows=nrow(sn),
Temp=(sn$Temp -19)/18 + 1,
Ammonia=sn$Ammonia,
pH=sn$pH,
Ecoli=sn$Ecoli,
Coliforms=sn$Coliforms,
Time=sn$Time,
OL=ifelse(sn$OL==.5, 1, ifelse(sn$OL==3.5, 3, sn$OL)),
Enterococci=sn$Enterococci,
sCOD=sn$sCOD,
VS=sn$VS,
TS=sn$TS,
Recipe=as.numeric(sn$Recipe),
Methane=sn$Methane
)
fit <- stan(
file = "./models/model_min.stan",
verbose = T,
chains=4,
data = stan_data
)
##
## TRANSLATING MODEL 'model_min' FROM Stan CODE TO C++ CODE NOW.
## successful in parsing the Stan model 'model_min'.
##
## CHECKING DATA AND PREPROCESSING FOR MODEL 'model_min' NOW.
##
## COMPILING MODEL 'model_min' NOW.
##
## STARTING SAMPLER FOR MODEL 'model_min' NOW.
## Warning: There were 1346 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
plot_pars <- names(fit)[!startsWith(prefix = "bugs", names(fit))]
plot(fit, pars=plot_pars, plotfun = "trace", inc_warmup = F)

mcmc_rhat(rhat(fit)[!startsWith(prefix = "bugs", names(fit))]) + yaxis_text(hjust = 0)

ggplot(data.frame(x = c(-5, 5)), aes(x)) +
stat_function(fun = dcauchy, n = 1e3, args = list(location = 1, scale = 1))

stan_dens(fit, pars=c("recipe_effect"))

plot(fit, pars=c("recipe_effect"))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

plot(fit, pars=c("temp_effect"))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

plot(fit, pars=c("lp__"))
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

#stan_model(fit)
(mle = optimizing(stan_model("./models/model_min.stan"), data=stan_data))
## $par
## recipe_effect[1] recipe_effect[2] recipe_effect[3] temp_effect[1]
## 1.786915e+00 -4.330820e-01 -6.778961e-01 4.311390e-02
## temp_effect[2] temp_effect[3] ol_effect[1] ol_effect[2]
## -9.153030e-01 1.673147e-01 7.821730e+00 1.626072e+01
## ol_effect[3] sigma bugs_hat[1] bugs_hat[2]
## 1.787536e+01 2.160753e+30 1.183533e+05 1.125921e+05
## bugs_hat[3] bugs_hat[4] bugs_hat[5] bugs_hat[6]
## 1.657377e+05 1.254627e+05 7.855682e+04 1.099506e+05
## bugs_hat[7] bugs_hat[8] bugs_hat[9] bugs_hat[10]
## 1.748567e+05 1.679054e+05 2.132778e+05 1.667928e+05
## bugs_hat[11] bugs_hat[12] bugs_hat[13] bugs_hat[14]
## 1.571769e+05 2.121249e+05 1.706268e+05 1.724419e+05
## bugs_hat[15] bugs_hat[16] bugs_hat[17] bugs_hat[18]
## 6.044487e+04 6.027080e+04 1.651412e+05 1.685152e+05
## bugs_hat[19] bugs_hat[20] bugs_hat[21] bugs_hat[22]
## 2.209633e+05 8.191879e+04 7.832521e+04 1.726854e+05
## bugs_hat[23] bugs_hat[24] bugs_hat[25] bugs_hat[26]
## 1.166730e+05 1.688998e+05 1.087502e+05 1.694282e+05
## bugs_hat[27] bugs_hat[28] bugs_hat[29] bugs_hat[30]
## 7.786855e+04 4.105228e+05 1.130727e+05 1.178734e+05
## bugs_hat[31] bugs_hat[32] bugs_hat[33] bugs_hat[34]
## 1.116322e+05 1.810452e+05 8.631004e+04 1.725937e+05
## bugs_hat[35] bugs_hat[36] bugs_hat[37] bugs_hat[38]
## 5.615964e+04 1.154729e+05 1.537116e+05 5.951931e+04
## bugs_hat[39] bugs_hat[40] bugs_hat[41] bugs_hat[42]
## 1.860855e+05 1.237718e+05 1.786968e+05 1.685509e+05
## bugs_hat[43] bugs_hat[44] bugs_hat[45] bugs_hat[46]
## 1.690371e+05 5.822546e+04 1.644123e+05 1.365916e+05
## bugs_hat[47] bugs_hat[48] bugs_hat[49] bugs_hat[50]
## 2.109678e+05 1.952756e+05 1.583445e+05 1.165913e+05
## bugs_hat[51] bugs_hat[52] bugs_hat[53] bugs_hat[54]
## 8.122379e+04 1.608862e+05 1.305319e+05 6.703936e+04
## bugs_hat[55] bugs_hat[56] bugs_hat[57] bugs_hat[58]
## 1.332220e+05 8.295257e+04 6.266193e+04 2.928694e+05
## bugs_hat[59] bugs_hat[60] bugs_hat[61] bugs_hat[62]
## 2.928694e+05 4.105228e+05 7.382467e+04 1.021236e+05
## bugs_hat[63] bugs_hat[64] bugs_hat[65] bugs_hat[66]
## 2.928694e+05 4.105228e+05 6.844350e+04 1.021236e+05
## bugs_hat[67] bugs_hat[68] bugs_hat[69] bugs_hat[70]
## 1.748564e+05 1.366878e+05 1.021236e+05 1.419258e+05
## bugs_hat[71] bugs_hat[72] bugs_hat[73] bugs_hat[74]
## 1.123880e+05 2.021011e+05 1.378997e+05 7.966518e+04
## bugs_hat[75] bugs_hat[76] bugs_hat[77] bugs_hat[78]
## 8.177432e+04 1.809906e+05 1.855631e+05 8.421338e+04
## bugs_hat[79] bugs_hat[80] bugs_hat[81] bugs_hat[82]
## 2.449770e+05 2.028970e+05 2.518740e+05 1.412786e+05
## bugs_hat[83] bugs_hat[84] bugs_hat[85] bugs_hat[86]
## 2.450736e+05 7.424722e+04 2.154774e+05 7.825904e+04
## bugs_hat[87] bugs_hat[88] bugs_hat[89] bugs_hat[90]
## 1.328324e+05 3.188351e+05 8.410832e+04 1.836829e+05
## bugs_hat[91] bugs_hat[92] bugs_hat[93] bugs_hat[94]
## 1.436913e+05 3.188351e+05 5.950899e+04 1.406585e+05
## bugs_hat[95] bugs_hat[96] bugs_hat[97] bugs_hat[98]
## 2.409977e+05 8.668223e+04 1.860511e+05 1.419259e+05
## bugs_hat[99] bugs_hat[100] bugs_hat[101] bugs_hat[102]
## 7.309299e+04 9.629382e+04 1.944402e+05 2.048183e+05
## bugs_hat[103] bugs_hat[104] bugs_hat[105] bugs_hat[106]
## 1.746205e+05 1.679254e+05 2.450736e+05 1.309019e+05
## bugs_hat[107] bugs_hat[108] bugs_hat[109] bugs_hat[110]
## 2.362202e+05 5.950849e+04 2.317106e+05 1.465055e+05
## bugs_hat[111] bugs_hat[112] bugs_hat[113] bugs_hat[114]
## 2.159638e+05 1.877536e+05 1.717681e+05 1.474174e+05
## bugs_hat[115] bugs_hat[116] bugs_hat[117] bugs_hat[118]
## 7.355702e+04 2.017443e+05 1.870239e+05 2.498922e+05
## bugs_hat[119] bugs_hat[120] bugs_hat[121] bugs_hat[122]
## 2.351592e+05 2.450736e+05 5.840009e+04 9.629382e+04
## bugs_hat[123] bugs_hat[124] bugs_hat[125] bugs_hat[126]
## 1.523874e+05 6.468478e+04 1.262270e+05 7.713658e+04
## bugs_hat[127] bugs_hat[128] bugs_hat[129] bugs_hat[130]
## 9.629382e+04 1.391061e+05 2.322411e+05 1.756472e+05
## bugs_hat[131] bugs_hat[132] bugs_hat[133] bugs_hat[134]
## 1.304193e+05 1.422432e+05 1.339400e+05 1.400714e+05
## bugs_hat[135] bugs_hat[136] bugs_hat[137] bugs_hat[138]
## 7.228446e+04 1.357281e+05 7.471261e+04 3.188351e+05
## bugs_hat[139] bugs_hat[140] bugs_hat[141] bugs_hat[142]
## 1.551771e+05 2.276386e+05 2.005904e+05 8.619428e+04
## bugs_hat[143] bugs_hat[144] bugs_hat[145] bugs_hat[146]
## 1.370592e+05 1.279139e+05 1.946281e+05 1.376996e+05
## bugs_hat[147] bugs_hat[148] bugs_hat[149] bugs_hat[150]
## 6.227995e+04 1.327276e+05 1.329684e+05 2.005906e+05
## bugs_hat[151] bugs_hat[152] bugs_hat[153] bugs_hat[154]
## 1.848310e+05 7.077057e+04 9.186715e+04 2.086660e+05
## bugs_hat[155] bugs_hat[156] bugs_hat[157] bugs_hat[158]
## 2.327458e+05 1.322461e+05 6.098742e+04 8.526667e+04
## bugs_hat[159] bugs_hat[160] bugs_hat[161] bugs_hat[162]
## 1.293572e+05 8.202186e+04 6.737795e+04 1.366524e+05
## bugs_hat[163] bugs_hat[164] bugs_hat[165] bugs_hat[166]
## 1.558638e+05 2.013722e+05 1.332089e+05 1.832947e+05
## bugs_hat[167] bugs_hat[168] bugs_hat[169] bugs_hat[170]
## 1.353303e+05 2.008834e+05 5.562575e+04 8.306325e+04
## bugs_hat[171] bugs_hat[172] bugs_hat[173] bugs_hat[174]
## 1.320042e+05 1.844476e+05 1.456484e+05 2.096003e+05
## bugs_hat[175] bugs_hat[176] bugs_hat[177] bugs_hat[178]
## 1.475041e+05 2.518917e+05 1.416325e+05 7.065319e+04
## bugs_hat[179] bugs_hat[180] bugs_hat[181] bugs_hat[182]
## 2.048193e+05 1.559704e+05 2.459603e+05 1.702286e+05
## bugs_hat[183] bugs_hat[184] bugs_hat[185] bugs_hat[186]
## 6.024708e+04 7.848485e+04 1.755901e+05 2.132718e+05
## bugs_hat[187] bugs_hat[188] bugs_hat[189] bugs_hat[190]
## 9.079141e+04 1.993639e+05 2.307861e+05 2.133430e+05
## bugs_hat[191] bugs_hat[192] bugs_hat[193] bugs_hat[194]
## 1.498855e+05 1.360102e+05 8.133175e+04 1.571327e+05
## bugs_hat[195] bugs_hat[196] bugs_hat[197] bugs_hat[198]
## 7.006891e+04 9.079141e+04 1.993639e+05 2.307861e+05
## bugs_hat[199] bugs_hat[200] bugs_hat[201] bugs_hat[202]
## 1.906687e+05 8.063218e+04 2.239904e+05 1.323661e+05
## bugs_hat[203] bugs_hat[204] bugs_hat[205] bugs_hat[206]
## 1.387446e+05 1.433174e+05 9.079141e+04 1.993639e+05
## bugs_hat[207]
## 2.307861e+05
##
## $value
## [1] -14603.89
##
## $return_code
## [1] 0